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SUMMARY 


A new technique for monitoring water levels in pressurized stainless steel cyl- 
inders has been developed. It is based on the differences in gamma ray attenuation 
coefficients in water and air. For an open-face, full-scale-model steel cylinder and 
a 10-p.Ci Cs ' gamma ray source, the counting rates for air and water differ by a 
factor of 15.60 ± 0.60 for cutoff energy of 200 keV and by a factor of 17.24 + 1.40 
for cutoff energy of 511 keV. These results are in reasonable agreement with the 
respective calculated values of 15.98 + 1.50 and 17.68 + 1.95. The respective values 
when the air is pressurized to 408.2 atm are calculated to be 4.75 ± 0.51 and 
5.01 ± 0.64. This large difference in counting rates ensures a clear capability for 
detecting the presence or absence of water in the gamma ray path in pressurized 
vessels. Computer programs for calculating nuclear radiation transmission are 
included in the appendix. 


INTRODUCTION 

Thermal probes in the Langley 8-Foot High-Temperature Tunnel have to be water- 
cooled in order to protect them from burnup during tunnel runs. The cooling water is 
stored in 30-ft-high, 16-in-diameter steel cylinders. Figure 1 shows a general view 
of the cylinder assembly for storage of cooling water. The water level in the cyl- 
inder is initially monitored using a simple U-tube manometer. Before the start of a 
tunnel run, high-pressure (about 408 atm) air is introduced into the cylinder above 
the water level and forces the cooling water through the multi jacketed thermal 
probes. Currently, there are no means available for directly monitoring the water 
level in the reservoir cylinders once the cylinders are pressurized. It is presumed 
that the initial manometer reading in the "nonpressurized" state gives correct indi- 
cation of the initial water level in the reservoir cylinder. This, of course, is 
true only if there was no residual pressurized air above the water level during the 
water pressure reading. It is therefore necessary to develop an independent tech- 
nique that monitors the presence of water at a critical height in the pressurized 
water reservoir cylinder and thus ensures an adequate supply of cooling water for the 
duration of the test. One technique that promises to meet this objective is based on 
differences in gamma ray attenuation coefficients in water and air. The principle of 
this technique as well as experimental results obtained with a full-scale mock-up 
model are described in the following sections. 


SYMBOLS 

c velocity of light 

D internal diameter of the steel reservoir cylinder 

E dis cutoff energy 

h Planck ' s constant 

I fc transmitted intensity 



incident intensity 


m Q rest mass of the electron 

r Q classical radius of the electron 

fraction by weight of the ith element 
x path length through a medium 

Z atomic number 

0 scattering angle 

|i total linear attenuation coefficient of the medium for the incident photons 

\i attenuation coefficient for the Compton process, given by the product of 

atomic density and Conpton scattering cross section per atom 

li linear attenuation coefficient for the ith element 

1 

ll attenuation coefficient for the photoelectric effect, given by the product 

r ne 

* of atomic density and the photoelectric cross section per atom 

\i attenuation coefficient for the pair production process, given by the 

^ product of atomic density and pair production cross section per atom 

v incident photon frequency 

0 

v' scattered photon frequency 

p density 

p density of the ith element 

1 

Q solid angle 

Subscripts: 
a air 

s steel 

w water 


PRINCIPLE OF OPERATION 

When a collimated beam of photons traverses a medium, it can interact with the 
atoms of the medium in a number of ways, depending upon its energy. These interac- 
tions result in the removal of the photon from the incident beam. The transmitted 
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intensity 1^ of a collimated beam of photons of incident intensity I after tra- 
versing through a thin absorber of thickness x may be written as follows (ref. 1): 


I 


t 


I e"^ 


( 1 ) 


where \i is the total linear attenuation coefficient of the medium for the incident 
photons. ji includes the effects of all types of interactions that the incident 
photons suffer in the absorber. For photons in the energy range 0.1 to 10.0 MeV, the 
dominant mechanisms are the photoelectric effect, the Compton effect, and the pair 
production effect; i.e.. 


M- = 


)i + |i + 
r pe r c 


PP 


(2) 


where 

|i = Attenuation coefficient for the photoelectric effect, given by the product of 
atomic density and photoelectric cross section per atom 

|i = Attenuation coefficient for the Compton process, given by the product of atomic 
density and Compton scattering cross section per atom 

|i = Attenuation coefficient for the pair production process, given by the product 
^ of atomic density and pair production cross section per atom 

Figure 2 shows the relative importance of these effects for various elements and for 
various photon energies (ref. 1). It is apparent that for Z < 40, the dominant 
mechanism for gamma ray energies of 0.2 to 7.0 MeV is the Compton scattering effect. 

The values of \x for different absorbers for photons in the energy range 0.01 
to 100.00 MeV have been reported by several authors (refs. 2 and 3). If the absorber 
is made up of a compound or a mixture of several elements, the overall attenuation 
coefficient is given by the following relation: 


if = 


P 




(3) 


where 

= Density of the ith element 

= Linear attenuation coefficient for the ith element 
w i = Fraction by weight of the ith element 
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If the absorber, on the other hand, is made up of a series of physically sep- 
arated materials, the transmission process through each material should be treated 
separately. Equation (1) in that case will take the following form: 
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( 4 ) 


where 

[ i ^ = Linear attenuation coefficient for the gamma rays in the ith component of 
distributed absorber 

= Path length in the ith component of the distributed absorber 

In the specialized case of a pressurized water vessel, equation (4) can be 
rewritten as follows: 


I (water) 


I e 


(5) 


I (air) = I e 
t o 




( 6 ) 


where 

|i = Linear attenuation coefficient for the gamma rays in steel 

LL = Linear attenuation coefficient for the gamma rays in water 
w 

\i = Linear attenuation coefficient for the gamma rays in air 

2x g = Gamma ray path length in steel (twice the steel cylinder wall thickness) 
3 ^ = Gamma ray path length in water (internal diameter of steel cylinder) 
x a = Gamma ray path length in air (internal diameter of steel cylinder) 
Combining equations (5) and (6), one obtains 


i t uir> 

I t (water) 


= e 


-( \x x -ji x ) 
a a w w 


I t (air) D ^ w -^ a ) 

■ ■ - — . = 0 
I fc (water) 
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where D is the internal diameter of the steel reservoir cylinder. The proposed 
technique is based on the fact that and ^ are very different for photons in 

the energy range 0.2 to 7.0 MeV. As indicated earlier, the linear attenuation coef- 
ficient of a medium is a very strong function of the photon energy. In order to 
obtain a reasonably high value of I (air) /I. (water) consistent with the reliability 

1 OT L. L -JOT 

of operation of the monitor, Cs was selected as the gamma ray source. Cs has a 
half-life of 30 years and emits a single gamma ray of energy 662 keV. As seen from 
figure 2, the dominant mode of interaction at this energy in all conponents of the 
pressurized stainless steel water vessel would be Compton effect. 

A theoretical discussion for the specific case of transmission of Cs 
(662 keV) gamma rays through the pressurized,* stainless steel, cooling water cylinder 
is given below. Figure 3 shows the essential geometry^jof the problem. 

Figure 4 illustrates the Compton scattering process for the case when the inci- 
dent photon energy far exceeds the binding energy of the struck electron. The dif- 
ferential Compton collision cross section per electron for incident unpolarized 
radiation of energy hv is given by the following equation (refs. 1, 4, and 5) 


da 

c 

d0“ 




sin 2 0l 


cm /electron 


( 8 ) 


where the scattered photon hv 1 goes into the solid angle dQ = 2rc sin 0 d9, and 
r Q is the classical electron radius (2.818 x 10~^ 3 cm). 

Because of the finite degree of collimation of the incident beam, finite thick- 
nesses of various components of the absorber medium in the path of the beam, and 
finite size of the radiation detector, multiple scattering effects must be taken into 
account (ref. 6). Because of the complexity of the secondary radiation produced in 
each collision, accurate theoretical calculation of broad beam attenuation is very 
difficult. Several different procedures have been used to calculate gamma ray trans- 
port through thick absorbers. Among the prominent methods reported (ref. 6) are 
(1) the method of successive scatterings, (2) the method of moments, and (3) the 
method of random sampling (Monte Carlo) . A Monte Carlo approach (refs. 6 and 7) was 
utilized in the present study in order to determine the number of photons, starting 
from a collimated Cs 137 source, that arrive at the uncollimated detector with an 
Energy > E dig . The theoretical details of the corrputational procedure are summarized 
in the appendix. The inclusion of multiple scattering effects reduces the difference 
in the counting rates between air and water media in the pressure vessel. For exam- 
ple, for air at a pressure of 1 atm, the ratio I t (air) /I t (water) falls from 20.17 
to 17.68 ± 1.95 when only scattered photons of Energy > 511 keV are counted and 
further falls to 15.98 ± 1.50 when the scattered photon energy discriminator level is 
reduced to 200 keV. For air at 408.2 atm, the corresponding values are 5.26, 

5.01 ± 0.64, and 4.75 ± 0.51, respectively. 


EXPERIMENTAL RESULTS 


A full-scale laboratory prototype system was constructed in order to test the 
concept of a nuclear technique for monitoring water levels. Figure 5 shows a sche- 
matic diagram of the experimental system for monitoring the presence or absence of 
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1 37 

water medium in the path of the test gamma rays. A 10-microcurie (nominal) Cs 
source provided the test gamma rays. The gamma rays transmitted through the pressure 
vessel were counted with a 2-in-diameter, 2-in-high Nal(Tl) crystal mounted on a 
high-gain photomultiplier. The photomultiplier output, after suitable amplification, 
was passed through a single-channel analyzer. Counts were recorded for predetermined 
intervals for two discriminator settings with and without water in the steel cyl- 
inder. The first discriminator setting of 200 keV was selected to provide for count- 
ing of all photons, having energies above the photoelectric energy limit, that 
arrived at the radiation detector. The counts registered will reflect considerable 
impact of multiple scattering effects in the absorbers between the source and the 
detector. The second discriminator setting of 511 keV was selected to restrict the 
counts to those photons that had suffered minimal scattering and, consequently, had 
lost minimal energy. The average values of I t (air) /I t ( water) for the two discrim- 
inator settings are summarized in table I. These values are in reasonable agreement 
with the corresponding calculated values listed in table II. 

Main reasons for the rather large errors in the experimental values lie in the 
weakness of the Cs 137 test source. For example, for E dis =511 keV, the counting 
rates with an empty cylinder with and without the source were approximately 
740 and 418 counts per minute (cpm) , respectively, whereas the counting rate with the 
source and with water in the cylinder was only 437 cpm. This poor signal-to- 
background ratio with water in the cylinder produces a large error in the 
I t (air) /I t ( water) value. A stronger test source should provide much better counting 
statistics and hence smaller error. This conclusion was verified by constructing a 
half-scale model and repeating all the aforementioned tests with it. The comparison 
between the experimental results and the calculated values is summarized in table 
III. The errors in the experimental results this time are considerably lower, as 
expected. 

Despite the large error in the I t (air ) /I t (water ) value with E dis = 511 keV 
in the full-scale model results, this discriminator setting would be preferable with 
a stronger radiation source because it makes the monitor more sensitive to the pres- 
ence or absence of water in the test radiation beam path. 

The prototype system, of course, could not be pressurized to 408 atm for reasons 
of safety. It was therefore not possible to determine experimentally the value of 
I t (air-pressurized)/l t (water) . However, there is no reason to doubt the theoret- 
ically calculated value, since the theory of gamma ray transmission through absorbers 
of various thicknesses is well understood. 

Figure 6 shows the experimentally observed counting rate for E^ g = 511 keV as 
the water was allowed to escape from the test cylinder. It is noted that the count- 
ing rate rises from an average value of 437 cpm to 740 cpm as the water is replaced 
by the air in the path of the test gamma rays. Allowing for a background counting 
rate of 418 cpm, these measurements give a value of approximately 17 for 
I t (air )/I t ( water) in reasonable agreement with the average value of 17.24 ± 1.40 
obtained under static conditions. 

In order to ensure that the tunnel could not be operated without adequate water 
supply in the pressurized water cylinders, a safety circuit based on the Cs^ 37 count- 
ing rate change has been designed, and a breadboard prototype has been built and 
successfully tested. The electronic safety circuit diagram is shown in figure 7. 

The circuit, besides providing a visible flashing light signal, will also provide an 
audio alarm signal if the water level is inadequate to start the tunnel operation. 

The pulses from the single-channel analyzer are shaped by the one-shot logic circuit 
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which increases the pulse width from 1 ps to 1.3 ms. This pulse width will permit 
handling counting rates as high as 769 counts per second (cps) . Higher counting 
rates will require shorter pulse width. The integration of these pulses produces a 
dc voltage proportional to the counting rate. The comparator circuit will turn on 
the alarm as this dc voltage matches or exceeds a preset reference level. 

CONCLUDING REMARKS 

The full-scale prototype system results under a pressure of 1 atm clearly demon- 
strate that the nuclear technique for monitoring water levels should easily distin- 
guish between water and air paths inside the water reservoir cylinders. Since the 
pressurized water cylinder tests could not be conducted with a laboratory prototype 
system, full-scale tests on actual cooling water cylinders were not performed. How- 
ever, reasonably good agreement between theoretical calculations and the experimental 
results for an open-faced prototype system indicates that equally good agreement 
should be obtained for a pressurized system, since the theory of electromagnetic 
radiation interaction with matter is well understood. 

The technique should, of course, be usable with other fluids, since the linear 
attenuation coefficients for intermediate energy gamma rays in air are considerably 
lower than in the fluids. It should also be adaptable to monitor continuously the 
fluid level in the reservoir systems and underground storage tanks. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 28, 1983 
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APPENDIX 


COMPUTER PROGRAMS FOR CALCULATING NUCLEAR RADIATION TRANSMISSION 

Symbols Used in Appendix 

c velocity of light 

D internal diameter of the steel reservoir cylinder 

d distance between the source and the detector 

ds path length between sampling planes 

E dis input cutoff energy 

Ef energy of the scattered photon 

Ej_ initial energy of the photon prior to scattering 

E q incident photon energy 

m Q rest mass of the electron 

n 0 product of the number of test atoms per cubic centimeter and the atomic 

number of the test atoms 

P c cumulative normalized probability 

P T total transmission probability 

P t transmission probability per photon 

R,p detector radius 

r Q classical radius of the electron 

Xf/y f ,z f x-, y-, z-coordinates following scattering 

x',y' ,z' final coordinates in the photon coordinate system 

2 

a = E /m c 

o o 

0 scattering angle 

9^ spherical coordinate defining gamma ray direction after scattering 

0^ spherical coordinate defining gamma ray direction before scattering 

9^ scattering angle in the photon coordinate system 

cMT) total scattering cross section 
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<!>£ spherical coordinate defining gamma ray direction after scattering 

spherical coordinate defining gamma ray direction before scattering 

4> scattering angle in the photon coordinate system 

s 

Q solid angle 


Program for Calculating Energy and Transmission 
Program Description 

The computer program RTRACK is written in FORTRAN IV language for the Control 
Data CYBER 170 series digital computer system with network operating system 
(NOS) 1.4. The program requires 31 600 octal locations of core storage. A typical 
case in which 5000 trajectories are computed requires approximately 350 CPU seconds 
on CYBER 173. 

The program models the stochastic process in which a gamma ray traverses a med- 
ium subject to multiple Compton scattering. When the detector location and size, 
initial gamma ray energy, dimensions and electron density of the intervening medium, 
and desired number of trajectories and path length for sampling are specified, the 
program uses a random number generator and an input probability distribution table to 
track the progress of each gamma ray. A separate program, described subsequently in 
this appendix, was developed to generate this table giving probability versus scat- 
tering angle ( 0 ) . 

Once the medium is defined, the z-axis is defined from the point source to the 
center of the detector with the origin at the source. Equally spaced planes are 
constructed perpendicular to this axis between the source and the detector with the 
spacing defined by the input sampling frequency parameter. The initial direction of 
the gamma ray is determined by a random value of <J>^ between 0 and 2n, and a random 
value of 9 ^ subject to collimation restrictions (where 4 ^ and 0 i are standard 
spherical coordinates) . An example of such a restriction is to limit the initial 
gamma ray to the forward hemisphere or to the cone defined by the source and detec- 
tor. A collimator with a diameter of 7/16 in. and a length of 1 27/32 in. was used. 

The gamma ray is allowed to travel along the initial trajectory until this tra- 
jectory intersects the first sampling plane. Hie transmission probability is com- 
puted using equation (A1) 


p t = 1 - a c (T > n e ds ( A1 ) 

where n 0 (given by the product of the number of test atoms per cubic centimeter and 
the atomic number of the test atoms, i.e. , n fi = nZ) is the electron density for the 
appropriate material (air, water, or the steel wall) , ds is the path length 
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between sampling planes, and a (T) is the total scattering cross section defined by 
equation (A2) 


0 .( 1 ) 



(A2) 


This transmission probability is compared with a uniformly distributed random number 
between 0 and 1 to determine if a collision occurs. If the random number is smaller, 
no collision occurs, and the gamma ray continues along the same path until it inter- 
sects the next sampling plane. Otherwise, a collision occurs, and a new path is 
determined. 

The scattering direction is determined by assigning a random value between 
0 and 2tz to <j) and by defining 0 g with a uniformly distributed random number 
and linear interpolation within the input probability distribution table. The 
z*-axis is defined as the direction of motion prior to the collision, and the new 
direction in the primed coordinate system is computed using the rotation matrices in 
equation (A3). 
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These coordinates are then transformed back to the laboratory coordinate system using 
the rotation matrices in equation (A4) . 
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where 9 and ^ define the gamma ray direction before scattering. The final 
direction in the laboratory coordinate system following the collision is then given 
by equations (A5) 



(A5) 


The gamma ray then travels along this trajectory until it intersects the next 
sampling plane. Following the scattering, the energy is reduced according to 
equation (A6) 





E - 

f 1 + a( 1 - cos 0 ) 


( A6 ) 


where is the initial energy and a = 

m c 2 
o 



This process is repeated for each sampling plane, and the gamma ray is said to 
reach the detector if, for the final plane, Xf 2 + y f 2 < R p 2 , where Rp is the 
detector radius. The gamma ray is dropped from further consideration if, for any 
sampling plane, one of the following inequalities is satisfied: x f > J d ' Yf > j d, 

z < 0, or E < E_. , where d is the distance between the source and the detector, 

and E^ is an input cutoff energy. 


After a gamma ray has reached the detector or has been dropped from further 
consideration, a new trajectory is initiated. This process is repeated until the 
desired number of trajectories has been completed. The total transmission probabil- 
ity is then defined by equation (A7) 


P 


T 


Gamma rays reaching the detector 
Total number of trajectories 


(A7) 


Description of FORTRAN Variables 

The following list contains a description of the significant FORTRAN variables 
appearing in the program. The dimension of each array is in parentheses beside the 
variable. Each variable is also identified as I, input variable; P, program vari- 
able; or 0, output variable. 


11 



APPENDIX 


FORTRAN Variable Type 

DIAM P 

DS I 

E FINAL I 

ENER P 

E2ER0 I 

NHIT 0 

NNE O 

NNX O 

NNY 0 

NNZ O 

NTRY I 

NZ I 

PHI P 

PHIF P 

PHIS P 

PTM P 

PTS P 

R P 

SEP I 

TDIAM I 

THETA P 

THETAF P 

THETAL P 

THETAS P 

THIK I 

TRANS O 


Description 

Diameter of detector, cm 

Spacing between sampling planes, cm 

Cutoff energy, keV 

Instantaneous gamma ray energy, keV 
Initial energy, keV 

Number of gamma rays counted by the detector 

Number of gamma rays with energy less than cutoff 
value 

Number of gamma rays exiting in x-direction 

Number of gamma rays exiting in y-direction 

Number of gamma rays exiting in z-direction 

Total number of trajectories 

Electron density of the interior medium, 
electrons/cm^ 

4> i of path before scattering 

<|)f of path after scattering 

<J) , scattering angle 
s 

Transmission probability for the interior medium 

Transmission probability for the steel wall 

Random number uniformly distributed between 0 and 1 

Distance between source and detector, in. 

Diameter of detector, in. 

0^ of path before scattering 

9 f of path after scattering 

Limiting value of 0^ due to collimation 

0 , scattering angle 
s 

Thickness of steel wall, in. 

Total transmission probability 
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TSCAT 
XFINAL 
XWALL1 
XWALL2 
XX( 180) 
YY ( 180) 


P Total scattering cross section 

P z-coordinate of the detector 

P z-coordinate of first steel/xnedium interface 

P z-coordinate of second steel/medium interface 

I Probability distribution table, 9 values 

I Probabililty distribution table, P Q values 


Program Listings 

Program RTRACK .- RTRACK, the main program, perforins all input and output opera- 
tions, controls the flow through the sampling planes, and calls the routines which 
perform the scattering transformations. A listing of RTRACK follows. 


1 


5 


10 


15 


20 


25 


30 


35 


40 


PROGRAM RTRACK ( OUTPUT , INPUT , TAPE 6-OUTPUT, TAP E5* INPUT , TAPE1 , 
1 T APE 2 ) 

COMMON PI, P0, ALPHA 
REAL NZDX,NZ 

DIMENSION XX(180),YY(180) 

DIMENSION ENN(5),C0N(5),PHN(5) 

EXTERNAL FUN 

DATA ENN/200.,300.,400.,500.,600./ 

DATA CON/10.57, 9. 19, 8 . 2 3, 7. 5 2, 6 . 96/ 

DATA PHN/ 2. 2 3, 0.66, 0.29, 0.16, 0.10/ 

1 F0RMAT(1H1,10X,5HRZER0,10X, 5HAL PHA, 9X, 6HVAL L 1,9X,6HWALL 2, 
1 7X,8HDISTANCE, 

1 4X, 11HTARGET D I AM, 2 X, 27HT0T SCAT CROSS SECTION (P T)///> 

2 FORMAT (1H1,5X,16HEL EC TR ON DENSITY, 5X, 15H PATH LENGTH , 

1 15X, 1HP/// ) 

3 F0RMAT(1H0,1PE15.4,1PE15.4,1PE15.4,1PE15.4,1PE15.4, 

1 1PE15.4,1PE29. 4 ) 

4 F0RMAT<//1H0,10X,*NZ * * , IP E 16 • 4 ) 

5 FORMAT (//1H0,14X,2HDS,8X, 6 HN Z DS ) 

6 F0RMAT(1H0,1PE16.4,1PE14.4) 

7 FORMAT ( *0 HITS, MISSES , X-0UT, Y-0UT, E-OUT, Z-0UT^ , 615 ) 

8 FORMAT (1H0,1PE21.4,1PE20.4,1PE16.4) 

PI • AC0$<-1.) 

THETAS - 0. 

PHIS « 0. 

THETAF - 0. 

PHIF - 0. 

C 

C READ PROBABILITY TABLE AND GENERATE NORMALIZED 

C CUMULATIVE DISTRIBUTION 

C 

DO 10 1-1,180 

RE AD ( 2 ) XX ( I ) , YY ( I ) 

10 CONTINUE 
SUM - 0. 

DO 20 1-1,180 
SUM - SUM + YY ( I ) 

YY ( I ) - SUM 

XX(I) - XX(I) ♦ PI/360. 

20 CONTINUE 

DO 30 1-1,180 
YY ( I ) • YY(I)/SUM 
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45 


50 


55 


60 


65 


70 


75 


80 


85 


90 


95 


100 


30 CONTINUE 

RO * 2. 818E-13 
40 READ(5>*) EZEPO, EFINAL 
IF ( EOF ( 5 ) .NE • 0 ) GO TO 999 
ALPHA - EZEPO/511. 

READ ( 5** ) THIK, SEP 
XINIT « 0. 

XVALL1 ■ 2.54*THIK 
XWALL1 - 2.+XWALL1 
XWALL2 * 2 • 54*( SEP - THIK) 

XVALL2 • 2.+XWALL 2 
XFINAL * 2.54*SEP 
XF INAL - 2 .*XFINAL 
PE AD ( 5> * ) TDIAM 
DIAM * 2 • 54*TDI AN 
TSCAT * GLEG15(0.,PI,FUN ) 

WR I TE ( 6, 1 ) 

WP ITE ( 1 , 2 ) 

WRITE (6,3) P0,ALPHA,XWALL1,XWALL2, XFINAL, DIAM,TSCAT 

RE AD ( 5 , * ) NZ, NTPY 

WR ITE ( 6, 6 ) NZ 

READ( 5,*) DS 

NZD X * NZ*D$ 

WRITE (6,5) 

WR ITE ( 6, 6 ) DS,N ZDX 
PTM • EXP(-TSCAT*NZDX) 

STELNZ * 22.07E+23 

PTS - EXP(-TSCAT*STELNZ*DS) 

NMISS - 0 
NHIT - 0 
NNX * 0 
NNY * 0 
NNE » 0 
NNZ ■ 0 
THETAL - PI/2. 

THETAL * ATAN (7 • /59. ) 

NSC AT * 0 
C 

C BEGIN HA IN LOOP OVER TRIALS 

C 

DO 150 1*1 , NTP Y 
ENER • F ZERO 
ZA * XINIT 
XA - 0. 

YA * 0. 

THETAI « RANF (T )*THETAL 
PHI I * RANF(T)*2.*PI 
THETA * THETAI 
PHI * PHI I 
50 CONTINUE 
R « RANF(T) 

PT * PTM 

I F ( Z A • L T . XWALL1 ) PT * PTS 

IF(ZA. GT.XWALL2.AND.ZA.LT. XFINAL) PT * PTS 
C 

C TEST TO SEE IF A SCATTER OCCURS 

C 

IF(R.LT.PT) GO TO 100 
C 

C YES, THERE IS A SCATTER 

C 
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105 


NSC AT * NSC AT + 1 
DO 60 K * 1 , 4 

IF(ENER.GE.ENN(K). AND. ENEP. LE.ENN(K^D) GO TO 70 
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110 


115 


120 


125 


130 


135 


140 


145 


150 


155 


160 


165 


60 CONTINUE 
K * 4 

70 CONTINUE 

CONX * CON ( K ) ♦ ( CON ( K+l ) - C ON C K ) ) ♦ ( ENEP - ENN (K ) ) 
1 /(ENN(K*1) - ENN (K ) ) 

PHNX - PHN(K) + ( PHN < K + 1 ) - PHN ( K ) >* < ENER - ENN < K ) ) 
1 / < ENN t K+l ) - ENN(K) ) 

FPACN « PHNX / (PHNX ♦ CONX) 

FPACN ■ 1. - FPACN 
R * 0# 

I F < ZA • LT • XWAL L 1 ) P « RANF(T) 

IF ( ZA *GT #XVALL2 ) P * RANF ( T ) 

C 

C TEST FOP PHOTOELECTRIC EFFECT 

C 

IF ( FRACN*LT*R ) GO TO 130 
C 

C COMPTON SCATTER 

C 

P - RANF(T) 

DO 00 K-1,179 

IF(YY(K).LE.R.AND*YY<K+1).GE.P) GO TO 90 
80 CONTINUE 
K « 179 
90 CONTINUE 

THFTAS « XX ( K ) ♦ (XX(K + 1) - XXCKMMR - YY ( K ) ) / 

1 ( YY ( K + 1 ) - YY ( K ) ) 

PHIS » PANF(T)*2^PI 
THETA « THETA*180./PI 
C 

C CONVERT ALL ANGLES TO DEGREES 

C 

PHI - PHI*180./PI 
THETAS * THETAS*180./PI 
PHIS - PHI S* 1 80 • / P I 
C 

C GET NEW DIRECTION 

C 

CALI POST < THETA, PHI, THE T AS , PHI S * THET AF , PHIF ) 

C 

C CHANGE ANGLES BACK TO RADIANS 

C 

THETA - THETAFtPI/180. 

PHI « PHIF*PI / 180. 

THETAS * THETAS*PI/1 80. 

C 

C DECREMENT ENERGY 

C 

ENER ■ ENEP / ( 1 • ♦ AIPHAM1* - COS ( THETA S )) ) 

100 CONTINUE 
C 

C INCREMENT PATH 

C 

ZA * ZA + DS*CO $ ( THETA ) 

XA * XA + DS*SIN(THETA)*COS (PHI ) 

YA « YA + DS*SIN(THETA)*SIN( PHI) 

C 

C TEST FOP OUTSIDE X-BOUNDS 

C 

IF(ABS<XA).GT.XFINALf 2. ) GO TO 110 
C 

C TEST FOP OUTSIDE Y-BOUNDS 
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170 


175 


180 


185 


190 


195 


200 


205 


210 

999 


C 

IF(ABS<YA).GT*XFINAL/2.) GO TO 120 
C 

C TEST FOR ENFRGY LESS THAN CUTOFF 

C 

IF(ENEP.LT.EFINAL) GO TO 130 
DZ * ABS ( XFINAl - ZA) 

C 

C TEST FOR OUTSIDE Z-BOUNDS 

C 

IF(ZA.LT*0.) GO TO 140 
C 

C TEST FOR PROXIMITY TO DETECTOR 

C 

IF ( DZ #GT • 4**DS ) GO TO 50 
DHIT * SORT < YA**2 ♦ XA**2> 

C 

C TEST FOP A HIT 

C 

IF(DHIT*LT*(DIAM/2# ) ) NHIT - NHIT + 1 
IF(DHIT.GE.(DIAM/2.)) NMISS * NMISS ♦ 1 
GO TO 150 
110 CONTINUE 

NNX • NNX ♦ 1 
NMISS * NMISS ♦ 1 
GO TO 150 
120 CONTINUE 

NN Y ■ NNY + 1 
NMISS « NMISS ♦ 1 
GO TO 150 
130 CONTINUE 

NNE * NNE + 1 
NMISS * NMISS ♦ 1 
GO TO 150 
140 CONTINUE 

NN Z » NNZ + 1 
NMISS « NMISS + 1 
150 CONTINUE 

WRITE (6*7) NHIT > NM IS S , NNX > NNY, NNE, NNZ 
TRANS « FLOAT (NTRY - NMISS )/ FLOAT ( NTRY ) 
WR I TE ( 1 * 8 ) NZfDSf TRANS 
GO TO 40 

CONTINUE 

STOP 

END 
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Function GLEG 15 .- Function GLEG 15 determines the value of a definite integral 
using a 15-point Gauss -Legendre quadrature scheme. It is used in determining the 
value of TSCAT. 


1 FUNCTION GLEG15 ( A, B > FUNCT ) 

C 

C FUNCTION GLEG15 EVALUATES INTEGRALS OF THE FORM 
C 

5 C B 

C INTEGRAL tG ( X )] OX 
C A 

C 

C WHERE A AND B ARE FINITE USING 15 POINT GA USS-L E GENDP E QUADRATURE. 
10 C 

C USE 


C 

C VALUE*GLEG15(A»P>FUNCT) 

C 

15 C A THE LOWER LIMIT OF INTEGRATION. MUST BE FINITE. 

C 

C B THE UPPER LIMIT OF INTEGRATION. MUST BE FINITE. 

C 

C FUNCT. .THE NAME OF AN EXTERNAL SUBROUTINE TO CALCULATE THE 

20 C VALUE OF G(X) AT SPECIFIED POINTS. 

C SUBROUTINE FUNCT ( F> N ) 

C F..CN INPUT THE N VALUES OF X 

C ON OUTPUT THE N VALUFS G(X) 

C THAT IS> Fm-GIFCI)] I« 1 ,...,N 

25 C N . .THE NUMBER OF FUNCTION EVALUATIONS. 

C THE NAME GIVEN FUNCT MUST APPEAR IN AN EXTERNAL 

C STATEMENT IN THE CALLING PROGRAM 

DIMENSION W(8)»R (8), F (1 5 ) 

DATA R/0. , .201194093997434, .394151347077563, 

30 + .570972172608538, .724417731360170* .848206583410427, 

* .937273392400705, .98799251 8020485/ 

DATA W/. 202578241925561, .198431485327111, .186161000015562, 

* .166269205816993, . 1 39570677926154 ,. 10 71 592204671 71 , 

* . 703660474 88 1 081 E-01, . 307532419961172 E-01 / 

35 C» ( B-A ) *. 5 

D- (B+A )*• 5 
F ( 1 ) *D 
DO 10 1-2,8 

s-r<i>*c 

40 F ( I ) -D+S 

10 F ( 1+7 ) -D— S 

CALL FUNCTt F, 1 5 ) 

D-W(1)*F(1) 

DP 70 1-2,8 

^5 20 D-D+W ( I )* ( F ( I ) + F ( 1+7 ) ) 

GLEG15«C*D 

RFTURN 

END 
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Subroutine FUN .- Subroutine FUN is required by GLEG 15 to compute the integrand 
of the definite integral, in this case the Klein-Nishina formula. 


1 SUBROUTINE FUN(F>N) 

DIMENSION F ( 1 ) 

COMMON PI, PO, ALPHA 
DO 10 I«1,N 

5 THETA » F C I > 

ANS « 2.*PI*P0**2*SIN(THETA) 

TMP1 ■ 1. ♦ ALPHA* ( 1 • - COS ( THFTA ) ) 
TMP2 « 1. + COS (THETA )**2 
TMP 3 • (1. - COS (THETA ) ) **2 
10 ANS ■ ANS*(l./TMP11**2*(TMP2/2. ) 

1 *(1. * ALPHA **2*TMP3/(TMP2*TMP1 ) ) 

F ( I ) « ANS 
10 CONTINUE 
RFTL'RN 

15 END 


Subroutine POST .- Subroutine POST computes the new direction of a trajectory 
following a scatter by defining and multiplying the appropriate rotation matrices. 


1 SUBPOUTINF P0ST(T1,P1,T2,P2,T3,P3) 

DIMENSION A ( A# A ) 

CALL UPOT 0C ( 1 >T2, A ) 

CALL UAPPLY(0.*0.»1.,A,X,Y,Z) 

5 CALL UROTOC ( 3, P 2, A ) 

CALL UAPPLY(X,Y»Z»A>X»Y»Z) 

CALL UP0T0C ( 1»T1» A ) 

CALL UAPPLY(X»Y»Z,A,X>Y»Z) 

CALL UP0T0C ( 3 > P 1 » A ) 

10 CALL UAPPLY(Y»Y»Z,A»X,Y»Z) 

P3- ATAN(X/Y) *180. /3.141 59265 
IF (X.LT.0..AND.Y.IT.0. ) P3-360.-P3 
IF (X.LT.O..AND.Y.GE.O. ) P3-180.-P3 
IF (X.GE.0..AND.Y.GE.0. ) P3-90.+P3 
15 IF (X.GE.0..ANP.Y.LT.0. ) P3 — P3 

R*SORT(X*X+Y*Y) 

IF (Y.LT.O.) R— R 
T3«ATAN(R/Z)*180./3. 14159265 
IF (R.LT.0..AND.Z.LT.0.) T3-180.-T3 
20 IF (R.LT.0..AND.Z.GT.0. ) T3—T3 

IF (R.GE.0..AND.Z.GE.0. ) T3-T3 

IF (R.GE.0..AND.Z.LT.0. ) T3-180.+T3 

RETURN 

END 
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Subroutines UAPPLY, UROTOC, UCLR, and UIDENT .- Subroutines UAPPLY, UROTOC, UCLR, 
and UIDENT are utility routines required by subroutine POST. 


1 SUBROUTINE UAPPLY (X,Y,Z» A, U,V,W) 

C 

C POST-MULTIPLIES A POINT VECTOR (X Y Z) BY 4*4 
C TRANSFORMATION MATRIX A TO GIVE (U V W). 

5 C 

REAL A(4,4) 

C H IS THE HOMOGENEOUS CO-ORDINATE » WHICH IS USED 
C TO EFFECT VARIOUS PROJECTIONS. IN THIS ROUTINE , 

C MATRIX A IS ASSUMED TO PRODUCE ONLY AFFINE 
10 C TRANSFORMATIONS. WHICH MEANS THE LAST COLUMN OF 

C MATRIX A IS (0 0 0 1), SO H IS UNITY. BUT IT,$ 

C HERE IN CASE ONE WISHES TO CHANGE IT... 

C (UNCOMMENT NEXT LINE AND DIVIDE X, Y, AND Z BY H.) 



C H « 

X*A (1*4) 

♦ Y*A ( 2, 4 ) 

♦ Z*A (3,4) 

♦ A ( 4, 4 ) 

15 

U ■ 

(y*A(l;l) 

♦ Y*A (2,1) 

♦ Z*A (3,1) 

♦ A (4,1) ) 


V * 

<X*A(1#2) 

+ Y*A (2,2) 

♦ Z*A (3,2) 

♦ A ( 4, 2 ) ) 


W - ( X*A ( 1,3) 

RFTURN 

END 

+ Y*A (2,3) 

+ Z*A (3,3) 

+ A ( 4, 3 ) ) 


1 


5 


10 


15 


20 


25 


SUBROUTINE UROTOC (IAXIS, ANGLE, A) 

C 

C CONSTRUCTS A 4*4 MATRIX THAT PERFORMS A ROTATION 
C ABOUT X-, Y— , OR Z-AXIS OF ANY ANGLE. 

C 


C INPUTS — 

C I AX IS 1«X 2-Y 3«Z 

C ANGLE AMOUNT OF ROTATION, IN CURRENT UNITS. 

C POSITIVE BY RIGHT-HAND PULE. 

C 

REAL A ( 4, 4 ) 

COMMON /TORADS/ ANGFAC 
DATA ANGFAC /. 01745329252/ 

C 


10 


IF ( I AXIS .LE. 0 .OR. IAXIS 
CALL UCLR (A) 

A ( 4 , 4 ) » 1. 

RAD - ANGFAC * ANGLE 
COSANG « COS (RAD) 

SINANG - SIN (PAD) 

GO TO (10, 20, 30), IAXIS 
A ( 1 , 1 ) « 1. 

A ( 2 , 2 ) » COSANG 
A ( 3, 2 ) ■ - SINANG 
A ( 2 , 3 ) * SINANG 
A ( 3 , 3 ) « COSANG 


.GE. 4) 


GO TO 40 


RETURN 

20 A ( 2 , 2 ) - 1. 

A ( 1 , 1 ) « COSANG 


MEASURED 
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30 


35 


AO 


30 


AO 


A ( 3 » 1 ) * SINANG 
A ( I » 3 ) * - SINANG 
A ( 3 » 3 ) * COSANG 
RETURN 
A ( 3 » 3 ) « 1. 

A C 1 » 1 ) » COSANG 
A(2»l) « - SINANG 
A(l>2) * SINANG 
A(2,2) * COSANG 
RETURN 

CALL UIDENT (A) 

RETURN 

END 


1 SUBROUTINE UCLR (A) 

C 

C CLEARS A* A MATRIX A TO ZEROS. 

C 

5 PEAL A ( A# A ) 

DO 10 1*1 » A 
DO 10 J*l» A 
10 A (I» J) - 0. 

RETURN 

10 END 


1 SUBROUTINE UIDENT (A) 

C 

C MAKES A+A IDENTITY MATRIX A. 
C 

5 REAL A(A,A) 

CALL UCLR (A) 

DO 10 I *1 # A 
10 A ( I » I ) - 1. 

RETURN 

10 END 
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Program Usage 

All program input is accomplished using FORTRAN list directed reads. The data 
may appear anywhere in the field, and when more than one item is specified, the data 
items are separated by commas. The first items in a data set are the initial gamma 
ray energy and the cutoff energy in keV. If, as a result of scattering, the gamma 
ray energy drops below the value of the cutoff energy, it is not counted by the 
detector. The next data items are the thickness of the container walls and the total 
distance between the source and the detector, in inches. These data are followed by 
the diameter of the detector in inches. These data items are illustrated in the 
sketch below: 


THICKNESS 


Cs 


137 


DETECTOR 

DIAMETER 


SEPARATION 


NoI(TI) 

DETECTOR 


The next data items include the electron density of the interior region (electrons/ 
cm 3 ) and the total number of gamma rays to be included in the run. The final data 
item is the path length to be used between sampling planes, in centimeters. If 
additional cases are desired in a single program execution, the entire data set must 
be repeated for each case. Sample data for a single program execution to analyze two 
cases where the interior region contains air and water are illustrated below: 


662. 0, 

200.0 


(EZERO, E FINAL) 

0.566, 

2.0 

8.0 


(THIK, SEP) 
(TDIAM) 

3.818E 

0.25 

+ 20, 

5000 

(NZ, NTRY) 
(DS) 

662.0/ 

0.566, 

2.0 

200.0 

8.0 


(Repeat for second case) 

3.346E 

0.25 

+ 23, 

5000 



In addition to the data described above, the program also requires a table con- 
taining probability versus scattering angle (9). These data should exist as an 
unformatted file on unit TAPE2. 
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The program generates formatted output on two files. File TAPE 6 contains a 
record of the program input as well as the number of gamma rays counted by the detec- 
tor. For those gamma rays not counted, TAPE 6 provides information as to whether they 
were lost because they exceeded the limiting dimensions in the x-, y-, or z- 
directions or they fell below the cutoff energy. File TAPE1 provides a brief summary 
of total transmission probabilities for each case considered. 


Program for Generating Scattering Angle Distribution 

Program TDIST2 was written to generate the probability distribution table 
required by program RTRACK. A listing of program TDIST2 is included in this section 
of the appendix. Note that function GLEG 15 and subroutine FUN, described with pro- 
gram RTRACK, are also required by this program. 

Program TDIST2 calculates the probability spectrum given by the following 
expression at intervals of 1°: 


da 

dQ 


1 


1 + a( 1 - cos 8) 


2 1 + cos 2 0) 


a 2 ( 1 - cos 0 ) 2 


( 1 + cos 9 ) [l + a ( 1 - cos 0 )] 


(A8) 


where r Q = 2.818 x 10 cm, and a = E 0 /m Q c^. This equation follows from equa- 
tion (8). Output from the program is written as an unformatted file to TAPE2 and is 
displayed in figure (A1). The normalized cumulative probability distribution 
versus 0 corresponding to this spectrum is shown in figure (A2) . Figure (A3) 
illustrates the results of using the distribution in figure (A2) to generate the 
frequency spectrum for 10 000 scattering events. As shown, both figures (Al) and 
(A3) are normalized to 10 000 events. 


1 


5 


10 


15 


20 


25 


PROGRAM TO I ST 2 ( OUTPUT *INPUT*TAPE6* OUTPUT *TAPE5* INPUT * 
1 TAPE2 ) 

COMMON PI*R0* ALPHA 
EXTERNAL FUN 
PI - AC0S (—1 . ) 

R0 * 2.818E-13 
ALPHA - 662 . /51 1 . 

TSCAT * GLEG15(0.*PI,FUN) 

WP ITE (6*1) 

1 FORMAT(1H1*10X*5HRZERO*10X*5HALPHA* 

1 2X* 27HT0T SCAT CROSS SFCTION (PT)///) 

WRITE (6*2) R0* ALPHA * TSC AT 

2 F0RMAT(1H0*1PE15.4*1PE22.4*1PE29.4) 

NUM « 180. 

00 100 1*1* NUM 

THETA « PI*FL0AT(I-1)/180. 

THETAP • THETA + PI/180. 

PSCAT * GLEG15(THETA*THETAP*FUN) 

WR ITE ( 2 ) THETA* PSCAT 
WP ITE ( 6* 3 ) THET A* PSCAT 

3 FORMAT ( 2E16. 8 ) 

100 CONTINUE 

CONTINUE 

STOP 

END 
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TABLE I.- SUMMARY OF EXPERIMENTAL VALUES FOR 
I t (air)/I t (water) FOR FULL-SCALE MODEL 


E dis' kev 

Air pressure, 
atm 

I fc (air) 

1^ (water) 

200 

1 

15.60 ± 0.60 

511 

1 

17.24 ± 1.40 


TABLE II.- SUMMARY OF CALCULATED VALUE S a FOR I fc ( air )/l fc ( water ) 
UNDER VARIOUS OPERATING CONDITIONS FOR FULL-SCALE MODEL 


(a) "Ideal" narrow beam geometry (no multiple scattering) 


E dis' keV 

Air pressure. 

V air) 

atm 

I (water) 

662 

1 

20.17 

662 

408.2 

5.26 


(b) "Poor" broad beam geometry (multiple scattering allowed) 


E dis' keV 

Air pressure, 
atm 

I fc (air) 

1^ (water) 

200 

1 

15.98 ± 1.50 

200 

408.2 

4.75 ± 0.51 

511 

1 

17.68 ± 1.95 

511 

408.2 

5.01 ± 0.64 


a These values were calculated for the following geomet- 
rical conditions: 


Collimator hole diameter = 0.438 in. 

Collimator depth = 1.844 in. 

Outer diameter of steel vessel = 16.000 in. 

Inner diameter of steel vessel = 13.736 in. 

Nal crystal dimensions = 2-in. diameter x 2-in. height 






TABLE III.- COMPARISON BETWEEN EXPERIMENTAL RESULTS AND CALCULATED VALUES FOR 

I t (air)/I t (water) FOR HALF-SCALE MODEL 



Air 

pressure, 

I t (air) 


E dis' kev 

1^ (water) 



CL Uill 

Experimental value 

Calculated value 

200 

511 

1 

1 

3.60 ± 0.02 
3.99 ± 0.02 

3.38 ± 0.23 
3.84 + 0.28 
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Figure 3.- Geometrical details of the gamma ray transmission computation problem. 

Linear dimensions in inches. 


(COMPTON ELECTRON) 



Figure 4.- Compton scattering effect. Paths of incident and scattered photons define 
scattering plane. Path of recoiling electron also lies in same plane. 
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COUNTS / 5 -MINUTE INTERVAL 



Figure 5.- Schematic diagram of experimental system for monitoring presence or 

absence of water in test cylinder. 


■WATER 


-WATER/AIR- 


(740 cpm) 


(437 cpm) 


I t (AIR) 

I t (WATER) 


- 17 


( Bkgd. = 418 cpm) 


RUNNING TIME SINCE START OF WATER DISCHARGE, MINUTES 


Figure 6.- Variation in counting rate as water is discharged from full-scale 
prototype water cylinder. (Ejj. =511 keV. ) 





















1. Report No. 2. Government Accession No. 

NASA TM-8565 1 

3. Recipient's Catalog No. 

4. Title and Subtitle 

DEVELOPMENT OF A NUCLEAR TECHNIQUE FOR MONITORING WATER 
LEVELS IN PRESSURIZED VESSELS 

5. Report Date 

September 1983 

6. Performing Organization Code 

505-33-53-14 

7. Author(s) 

Jag J. Singh, William T. Davis, and Gerald H. Mall 

8. Performing Organization Report No. 

L- 15642 

10. Work Unit No. 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, VA 23665 

11. Contract or Grant No. 

13. Type of Report and Period Covered 
Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 

14. Sponsoring Agency Code 

15. Supplementary Notes 


Jag J. Singh and William T. Davis: Iangley Research Center, Hampton, Virginia. 

Gerald H. Mall: Corrputer Sciences Corporation, Hanpton, Virginia. 


16. Abstract 


A new technique for monitoring water levels in pressurized stainless steel cylinders 
has been developed. It is based on differences in attenuation coefficients of water 
and air for Cs (662 keV) gamma rays. Experimentally observed gamma ray counting 
rates with and without water in model reservoir cylinders have been compared with 
corresponding calculated values for two different gamma ray detection threshold ener- 
gies. Calculated values include the effects of multiple scattering and attendant 
gamma ray energy reductions. The agreement between the measured and calculated val- 
ues is reasonably good. Computer programs for calculating angular and spectral dis- 
tributions of scattered radiation in various media are included. 


17. Key Words (Suggested by Author(s)) 18. Distribution Statement 


Gamma rays Unclassified - Unlimited 

Narrow beam geometry 
Broad beam geometry 
Attenuation coefficient 

Photoelectric effect Subject Category 35 

Compton scattering cross section 

Pair production cross section 

Multiple scattering 

Radiation detectors 

Monte Carlo simulation 

Pressurized vessels 

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of Pages 22. Price 

Unclassified Unclassified 33 A03 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1983 





















National Aeronautics and 
Space Administration 

Washington, D.C. 

20546 


Official Business 

Penalty for Private Use, $300 


THIRD-CLASS BULK RATE 


Postage and Fees Paid 
National Aeronautics and 
Space Administration 
NASA-451 





POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 



